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PREDICTABILITY OF THE EARTH’S POLAR MOTION 


B. Fong Chao 
ABSTRACT 

The numerical prediction of the Earth’s polar motion is of both theoretical and practical 
interests. The present paper is aimed at a comprehensive, experimental study of the predictability 
of the polar motion using a homogeneous BIH (Bureau International de l’Heure) data set for the 
period 1 967-1983. Based on our knowledge of the physics of the annual and the Chandler wobbles, 
we build the numerical model for the polar motion by allowing the wobble periods to vary. Using 
an optimum base length of 6 years for prediction, this “floating-period” model, equipped with a 
non-linear least-squares estimator, is found to yield polar motion predictions accurate to within 
0".012 to 0".024 depending on the prediction length up to one year, corresponding to a predict- 
ability of 91-83%. This represents a considerable improvement over the conventional fixed-period 
predictor, which, by its nature, does not respond to variations in the apparent wobble periods (in 
particular, a dramatic decrease in the periods of both the annual and the Chandler wobbles after 
the year 1980). The superiority of the floating-period predictor to other predictors based on 
critically different numerical models is also demonstrated. 
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PREDICTABILITY OF THE EARTH’S POLAR MOTION 
L. INTRODUCTION 

The art of prediction of a time series consists in the extrapolation into the future of past 
behavior of the dynamic system that generates the time series. Some systems are deterministic and 
hence completely predictable (at least in principle), such as tides or the time of the next sunrise. 
Other systems are for the most part statistical in nature and often too erratic to be subjected to 
accurate predictions; examples include the weather pattern and stock market tradings. 

The Earth’s polar motion is a dynamic system that lies somewhere between the two extremes 
(a quantitative treatment is given in Section 4). To first approximation, the polar motion is com- 
posed of a secular motion and two quasi-circular components: the annual wobble and the 14-month 
Chandler wobble. These components are by no means completely deterministic (as are, say, the 
tides), yet each of them, governed by different physics, is coherent enough with respect to time 
that a fairly accurate prediction can be expected. This makes the problem of polar motion predic- 
tion meaningful and, in my opinion, challenging. 

Besides theoretical interests, the prediction of polar motion is of practical value as well, as 
attested by the independent publication of predicted pole positions by two time services: the 
BIH (Bureau International de l’Heure) Circular B/C [Feissel, 1980] , and the USNO (U.S. Naval 
Observatory) Series 7 [Babcock, 1982] . The conventional prediction method adopted is one that 
can be characterized by fixed wobble periods. Zhu [ 1982] has made a comprehensive study of 
polar motion prediction that uses such a “fixed -period” scheme; it, as far as I know, has remained 
the only published work on the subject to date. 

The present paper is an experimental study of the predictability of the polar motion based on 
a homogeneous BIH data set for the period 1967-1983. It will be shown that with a floating- 
period numerical model and a capable estimation method, predictions accurate to within 0".01 2 
to 0".024 can be achieved for the polar motion, depending on the prediction length up to one year. 
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The superiority of this “floating-period predictor” to other predictors that are based on critically 
different numerical models (including the conventional fixed-period predictor) will be demon- 
strated. 

2. PREDICTION (using the Floating-Period Predictor) 

The procedure of numerical prediction consists of three stages: model identification, param- 
eter estimation, and extrapolation [Box & Jenkins, 1970] . 

2.1 Model Identification 

Based on our knowledge of the physics of the polar motion, we decompose the latter into 
three major components: 

polar motion = (secular motion) + (annual wobble) + (Chandler wobble) (1) 

For the purpose of prediction, we shall consider record segments with lengths shorter than, say, 8 
years. It is with respect to this time scale that we discuss the numerical representation of each 
component in equation (1). 

The secular motion (relative to the Conventional International Origin, or CIO) can be well- 
represented by a linear term. It is presumably the combination of a true polar wander and any 
periodic component with a period much longer than, say, 8 years. 

The annual wobble is believed to be driven largely by the variations in the atmospheric mass 
distribution [see e.g. Lambeck, 1980] . The solar energy received by the atmosphere cycles every 
365.25 days with a nearly constant peak-to-peak amplitude. However, the Earth’s meteorological 
transfer function being intricate and non-linear, it is only natural that the annual wobble thereby 
generated will have an amplitude that varies from year to year and a period that fluctuates about 
the driving period 365.25 days. 

The Chandler wobble is the Earth’s free oscillation that corresponds to the Eulerian nutation 
of a rigid body; but the exact physical nature of its variations with time is largely unknown. One 
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hypothesis statesthat the Chandler wobble is simply a single-period oscillation convolved with a 
sequence of random or non-random excitations. An alternative hypothesis calls for multiple (and 
fixed) natural periods [see.e.g. Colombo & Shapiro, 1968] . Both hypotheses give rise to (apparent) 
instantaneous periods and amplitudes that are slowly varying with respect to time, as manifested by 
actual observations. (For the single-component hypothesis, a good synthetic example is given by 
the “disturbed pendulum” model of Yule [ 1 927] ; for the modulating effects of multiple periods on 
the instantaneous period and amplitude of the polar motion, see Chao [1983] .) 

The above discussions justify the following numerical model for the polar motion: 

X(t) = A + Bt + C a cos^~+ + C c cos (p+ $cj (2) 

for the X-component, and a similar expression for the Y-component Y(t). A+Bt is the linear secular 
term, and the two sinusoidal terms represent the annual and the Chandler wobbles (C stands for the 
amplitude, P the period, and <p the phase angle). The eight (real) parameters A, B, C a , P a , 0 a , C c , 

P c , <f> c are unknowns to be estimated from a given record segment X(t) or Y(t) — each record seg- 
ment yields one set of estimated parameters. Since we allow the data in each given record segment 
to “choose” the wobble periods, we call equation (2) the floating-period model. Note that here 
we treat X(t) and Y(t) as independent series so as not to restrict ourselves only to the prograde 
portion of the true polar motion, as would if we chose to model the complex series X(t)+iY(t) in 
terms of “complex wobbles” of the form exp(i27rt/P+0). Although the physics behind equation 
(2) is for the most part unknown, we shall show that this model is adequate for the purpose of 
prediction. 

2.2 Parameter Estimation Method 

The next stage is, of course, the numerical estimation of the eight free parameters in equation 
(2). This is a non-linear fitting problem, for equation (2) is non-linear in some of the parameters 
(particularly the periods). Thus we implemented the non-linear least-squares estimator, called 
“CURFIT”, according to Bevington [ 1969, pp. 23 2-242] . It uses the Marquardt algorithm which 
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combines a gradient search with an analytical solution developed from linearizing the fitting func- 
tion. This algorithm is found to be robust and fast. 

2.3 Extrapolation and Error Analysis 

As an experiment, we now apply the floating-period estimator described above to the X- and 
Y-series of a homogeneous BIH polar motion data set, given at 5-day intervals spanning 1967.0- 

1983.3 with the origin adjusted to the CIO in 1967 [Feissel, 1980] . These two time series, hence- 
forth called BIH-X(t) and BIH-Y(t) respectively, are referred to the 1980 Nutation Theory [Capi- 
taine & Feissel, 1983] , and have been smoothed by means of the Vondrak algorithm. They are 
displayed in Figure 1 . 

Our basic procedure is straightforward. We simply “chop” an N-year segment (where N=5, 
6,7,8) from BIH-X (or BIH-Y), apply the floating-period estimator, and use the derived model 
parameters to extrapolate the time series M days (where M=25,50,100,180,365) into the future. 

We start from the first N-year segment at the beginning of BIH-X (or BIH-Y), make an M-day 
prediction, and then shift forward M days for the next M-day prediction. This procedure con- 
tinues across the length of BIH-X (or BIH-Y) until we hit the end. Each M-day prediction can now 
be compared with actual observations, and the root-mean-square (rms) prediction error (for each 
M-day prediction), denoted by E, can be evaluated. Finally, we can average the sequence of the 
prediction errors E obtained consecutively across the whole length of BIH-X (or BIH-Y) to yield 
the average prediction error E. The value of E will be used as a performance indicator for the pre- 
dictor under consideration. The prediction scheme described above will be denoted as [N(years), 
M(days)] , where N is referred to as the base length and M the prediction length; and we shall be 
speaking of E[N,M] (a sequence) and of E[N,M] (a single quantity). 

Figure 2 summarizes the performance of our floating-period predictor by showing E[N,M] for 
N=5~8(years) and M=25~365(days). Two things are to be noted: 

(i) Minimum E invariably occurs around N=6(years). This is in fact expected, because 6 years 
is the beating period between the annual and the Chandler wobbles. For N<6, the presence of 
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noise in the data greatly decreases the capability of our floating-period estimator (or any other 
estimator for that matter) in resolving the two wobbles. (Note, incidentally, that all we need 
to estimate the set of 8 parameters are 8 data points, which in the present case amount to a 
base length of only 40 days, if the data were noise-free.) For N>6, on the other hand, E in- 
creases because the coherency of each component in the polar motion decreases as time span 
gets large; and Figure 2 clearly implies that the coherence length of the polar motion is less 
than 6 years. Thus, we conclude that 6 years is the optimum base length for the prediction 
of polar motions. 




BASE LENGTH N (years) 


Figure 2. The average prediction error E[N,M] by the floating-period predictor, where N is the 
base length (in years) and M is the prediction length (in days). 


6 




(ii) For any base length N, E* increases monotonously as the prediction length M gets longer. 
This is again to be expected because of the decaying of the coherence, and hence predicta- 
bility, of the polar motion with respect to time. 

As an example, we plot in Figure 3 a typical 3-year segment (1977.7-1980.7) of observed pole 
positions (solid line) compared with those predicted by some ten consecutive [6,100] floating- 
period predictions (dotted line, at 5-day intervals). The two sequences E[6,100] (for BIH-X and 
BIH-Y) are shown as the dashed lines in Figure 4 (their averages, E[6,100] , are shown in Figure 
2 as well as in Table 1). Table 1 lists"E[6,M] vs. M, and is plotted in Figure 5 as dashed lines. From 
Table 1 we see that, depending on the prediction length (within 25~365 days), the floating-period 
predictor in general achieves accuracy to within 0".01 2 to 0".024. In comparison, the 6-year rms 
fit residuals of the floating-period estimator are in the range 0".01 1-0". 01 4. 



X (arcsec) 

Figure 3. A typical 3-year segment of the polar motion. The solid line is the BIH pole path, the 
dotted line shows the pole path predicted by the floating-period [6,100] scheme. The 
cross + gives the position of the 1967 CIO. 
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Table 1 

The average prediction error E[6,M] (in unit of 0".001) by the floating-period predictor, vs. the 

prediction length M (see dashed lines in Fig 5). 


M (days) 

25 

50 

100 

180 

365 

BIH-X 

15.1 

16.9 

19.6 

22.6 

23.6 

BIH-Y 

12.8 

13.9 

15.8 

17.5 

19.5 


3_. OTHER PREDICTORS 

In this section we shall compare other predictors which, when subjected to the same experi- 
ment described in Section 2.3, give less satisfactory performance than the floating-period predictor. 

3.1 Fixed-Period predictor 

Fixed-period predictor is one that uses pre-determined annual and Chandler wobble periods 
in the numerical model, hence the name “fixed-period”. Let the periods be 365.25 and 435.00 
days, respectively: 

/ 27rt \ / 2irt 

X(„ - A ♦ B. ♦ C,co^— ♦ *,) ♦ C c cos^— ♦ * 

and similarly for Y(t). This is the simplest polar motion model, and the estimation of the six 
parameters is a linear problem (the non-linearity with respect to 0 in expression (3) dissolves when 
we expand the cosine terms). To date, the fixed-period predictor has been the only predictor in 
use [BIH Rapid Service, USNO Series 7] , or discussed [Zhu, 1982] . 

However, when we apply the fixed-period predictor to the BIH data set, something unexpected 

appears. The solid lines in Figure 4 show the prediction error sequence E [ 6,1 00 ] by the fixed- 

period predictor for BIH-X and BIH-Y. Prior to the year 1980, this predictor worked steadily, 

although not quite as well as the floating-period predictor. For example, E[6,100] prior to 1 980 is 

found to be 0".021 for BIH-X and 0".018 for BIH-Y. This agrees with the results of Zhu [ 1982] , 

and is about 10% higher than those of the floating-period predictor (c.f. Table 1). However, 

E[6,100] from the fixed-period predictor soars to an average level of 0”.070 for BIH-X and 0 /, .059 

for BIH-Y after 1980 (see solid lines in Figure 4). The same phenomenon occurs in all other pre- 
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diction schemes [6,M] with different prediction lengths M. The root of this phenomenon is traced 
to a dramatic decrease in both annual and Chandler wobble periods starting around 1977 (see 
Section 4), a situation that the fixed-period predictor cannot cope. This is also reflected in a 
similar post-1980 jump in the 6-year fit residuals by the fixed-period estimator. 



Figure4. Comparison of the prediction error sequence E [ 6,1 00] by the fixed-period predictor 
(solid lines) and by the floating-period predictor (dashed lines). 
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The U.S. Naval Observatory (USNO), in its weekly “Time Service Announcement Series 7 
(Earth Orientation Bulletin)”, publishes 40 days in advance pole positions predicted by a fixed- 
period predictor. The major difference between the USNO scheme and the fixed-period scheme 
discussed above is in the base length used in prediction — USNO uses a base length of only 2 years 
[McCarthy, personal communication, 1983] . The advantage of using a short base length is that the 
resultant estimates are more “instantaneous”. In fact, our experiments have confirmed USNO 
Series 7 by showing that the fixed-period E[ 2,40] lies as low as 0".015, which is comparable to 
E[6,40] of our floating-period predictor (see the dashed lines in Figure 5). However, this success 
is only apparent. Indeed, in the fixed-period [2,M] scheme, since the base length falls much shorter 
than 6 years, the estimates are prone to noise in the data (see Section 2.3). In particular, this is 
found to result in wobble amplitude estimates that are unstable and sometimes unreasonable. For 
short prediction lengths (M shorter than 40 days, say) this creates little discomfort as far as the 
prediction is concerned because of the general compensating capability of least-squares estimators. 
However, when one extrapolates these (unstable) estimates further into the future, one would 
expect a rapidly growing deviation of the predictions from actual observations. This is indeed 
found to be the case - for example, E] 2,365] is found to be 0''.029, almost double the value of 
E[2,40] . Furthermore, although not as severely as discussed earlier, the fixed-period [2,M] scheme 
also sees a post-1980 jump in E[2,M] as M increases. 

In conclusion, the fixed-period predictor yields parameters that fail to reflect true variations 
in the polar motion; and, as a result, performs less satisfactorily than the floating -period predictor 
and is potentially dangerous in the sense that it may sometimes lead to disastrous predictions. 

3.2 Floating-Complex Period Predictor 

The superiority of the floating-period predictor to the fixed-period predictor is obviously the 
consequence of two additional degrees of freedom in the former fitting procedure, namely the two 
wobble periods. Similarly, we can build in yet more degrees of freedom to allow for variations in 
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the amplitude for each N-y ear segment. A natural choice is to make the wobble periods complex 
in the numerical model, so that their imaginary parts give exponential decay or growth in the 
wobble amplitudes, hence the name “floating-complex period”: 


/27rt \ /2nt 

X(t) = A+Bt + C a exp(a a t) cos(— + 0 a ) + C c exp(a c t) cos(-^— + 4> 

P a 


Vc 


( 4 ) 


and similarly for Y(t). Now we have ten (real) parameters to be estimated, some of which are non- 
linear; and again we use the “CURFIT” algorithm as in Section 2.2 


As expected, equation (4) is indeed found to be a better fitting model. For instance, the rms 
fit residual for each 6-year record segment.is invariably lower (in some cases as much as 30%) than 
that of the floating-period estimator. As a predictor, like the floating-period predictor, it also 
responds to the post-1977 drop in the wobble periods mentioned in Section 3.1 and shows no 
increase in the prediction errors E[6,M] after 1980 as did the fixed-period predictor. However, the 
average prediction error E[6,M] (M=25~365) of its predicted pole positions (by means of extra- 
polation) are generally higher than those by the floating-period predictor, as shown in Figure 5. 

This, of course, is indicative of the fact that the exponential decay/growth trend in the wobble 
amplitudes in general does not continue past, say, 6 years. Indeed, physically, we should not expect 
otherwise because the variations of wobble amplitude with respect to time are determined by the 
excitation functions which presumably are either erractic or changing from year to year. In con- 
clusion, we state that while the floating-complex period model does a good fitting job and may 
hence have more potential usages, it is less satisfactory as a predictor for the polar motion. 
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PREDICTION LENGTH M (days) 


Figure 5. Comparison of the average prediction error E[6,M] (as a function of the prediction 
length M) by the floating-complex period predictor (solid lines) and by the floating- 
period predictor (dashed lines, see Table I ). 


3.3 Autoregressive (AR) Predictor 

A related prediction scheme (but under a completely different philosophy) is the so-called 
linear prediction, which has been widely pursued in research fields such as geophysics, speech proc- 
essing, control theory, and economics [see e.g. Makhoul, 1975]. Its application in time-series 
prediction has been thoroughly studied by Box & Jenkins [1970] . Here we shall briefly develop 
one widely used class of numerical models in the linear prediction theory, namely the autoregres- 
sive (AR) models, and relate it to the present problem of polar motion prediction. 


An mth order AR model of a time series x(n) is 


x(n) = Sj x(n-i) + a(n) 

i= 1 

where] a(n)} is a random (or white) series, andjSp i=l ,2,. . .,m | , called the AR coefficients, are 


12 



the (real) model parameters yet to be determined. While having the physical meaning of being the 
external disturbances to the dynamic system, the series a(n) is just the nth prediction residual as far 
as prediction is concerned [for details see Box & Jenkins, 1970] . Many physical systems can be 
modeled by such an AR process; and the model identification problem consists mainly in the 
determination of the order m. Although the latter problem may often be different or even con- 
troversial, we shall show presently that it poses no obstacle to the modeling of polar motion time 
series. It is well-known that a pure sinusoidal function can be represented as a second-order AR 
time series [see e.g. Yule, 1927] . In fact, it can be shown that a linear combination of K complex 
sinusoids (with exponentially decaying or growing amplitudes) can be expressed as a 2Kth order AR 
time series, and the K complex periods are related to the 2K (real) AR coefficients through the 
Prony’s relation [for details see Chao & Gilbert, 1980; Chao, 1984]. For example, with a(n) 
representing some random excitation function and m=2, equation (5) is a mathematically simple, 
and physically sound model for the Chandler wobble. For a linear combination of two wobbles 
(K=2), namely the annual and the Chandler wobbles in the polar, motion, we let m=4. Equation 
(5) is then a general AR model which, apart from the secular term A+Bt, encompasses equation 
(2), (3), and (4). 

Now that the numerical model has been built, the estimation of parameters (i.e., the AR 
coefficients | Sj ;i= 1,2 ,3 ,4 1 ) is just a linear least-squares procedure, and the extrapolation is 
straightforward, all according to equation (5) by letting a(n)=0. In practice, unfortunately, we 
encounter the following two difficulties: 

(i) Since equation (5) does not allow for a linear term A+Bt, the latter has to be removed from 
the data prior to the AR estimation. This can be done by a least-squares fit. However, for 
record segments as short as several years, the fit is biased in the presence of the wobble terms. 
Synthetic experiments have indicated that this, in turn, can lead to biased estimations for the 
wobble parameters. (Nevertheless, it should be pointed out that the AR estimator is accurate 
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and powerful in analyzing zero-mean, trendless polar motion data, as demonstrated by Chao 
[1983]). 

(ii) A more serious problem which is inherent in the AR predictor results from the noise in 
the data. Thus, even with an accurate set of model parameters, a predicted value can be, and 
often is, seriously biased by the noise in the existing values used in the prediction. Put differ- 
ently, the noise in x(n— i) on the right side of equation (5) “propagates out” to x(n) on the left 
side of (5). Synthetic experiments show that average prediction errors E generated by an AR 
predictor are generally several times larger than the corresponding E’s by a floating-period 
predictor. 

Thus, although it is an adequate representation for the wobbles, the AR model (5), and hence the 

AR predictor, is unworthy in the prediction of polar motions. 


A CONCLUSIONS AND DISCUSSIONS 

As Zhu [ 1982] pointed out, the prediction errors consist of modeling errors and observation 
errors in the data. Assume that the modeling error and the observation error in both X- and Y- 
components are uncorrelated with one another. Then the sum of their individual variances equals 
the variance of the total prediction error: 


0^—0^ + o ^ 

°E .M.E. u O.E. 


( 6 ) 


where each term is in tum the sum of the variances for X- and Y-components ( a denotes the stand- 
ard deviation). Now let us define the predictability P simply as 


P = 1 - a 


M.E. 


I O 


P.M. 


(7) 


where o pjv1 is the standard deviation of the polar motion. For a system being modeled reason- 
ably we expect P to lie between 0 and 1 ; P=1 signifies a completely deterministic system which 
has been modeled exactly, while P=0 results from a complete failure in prediction with the 
modeling error as great as the signal itself. Obviously P depends on the predictor and the prediction 
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scheme in use, and represents a (linear) performance indicator. The quantity o p M can be estimated 
by computing the rms amplitude of the polar motion time series, and is found to be 0". 1 7. The 
quantity a M E can be determined using equation (6), where (i), o Q E is estimated to be about 
0".01 1 for our smoothed BIH data set [Feissel, personal communication, 1983] , and (ii) a £ can be 
estimated simply using E”, an essemble average of E which is itself an estimate (as a temporal mean) 
of a £ . Thus, from Table 1 and equation (7), we conclude that the predictability of the polar 

f 

motion using the floating-period predictor for prediction lengths in the range 25~365 days is 
91-83%. The predictability of the polar motion using the floating-complex period predictor is 
somewhat lower (90-77%); and those using the fixed-period and AR predictors are considerably 
lower yet, and, in the former case, sometimes ill-behaved (see Section 3.1 ). 

One obvious by-product of our prediction procedure is the set of sequences of estimated polar 
motion parameters as functions of time. One such sequence which is interesting and relevant to our 
present study is that of the estimated wobble periods. Figure 6 shows the variations in the esti- 
mated annual and Chandler wobble periods resulted from the [6,100] prediction scheme. The 
points represent 6-year running estimates of the periods, each shifted 100 days down the length of 
the BIH data set; thus features with scale shorter than, say, 2 years might be spurious (being simply 
the results of noise in the data). However, the sudden decrease starting around 1977 in both annual 
and Chandler wobble periods, which show up in both BIH-X (solid lines) and BIH-Y (dashed lines), 
is believed to be real. Other estimation techniques, namely the floating-complex period and the AR 
estimators, also unequivocally show a very similar trend in the period variations. This has caused 
our implementation of the fixed-period predictor to “fail” after 1980 as we have seen in Section 
3.1 . Presently we know nothing about its cause(s) or its geophysical implications. We can only 
speculate that this phenomenon has its root in the atmospheric circulation system; but this awaits 
further investigations. 
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Figure 6. Six-year running estimates of wobble periods estimated by the floating-period estimator 
from BIH-X (solid lines) and BIH-Y (dashed lines) for (a) the Chandler wobble, and (b) 
the annual wobble. The curves have been connected smoothly. 
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the Chandler wobbles after the year 1980). The superiority of the floating-period predictor to 
other predictors based on critically different numerical models is also demonstrated. 
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